Mon. Not. R. Astron. Soc. OOP. [TlfTSl (???) Printed 13 October 2011 (MN WT^ style file v2.2) 



On Point Spread Function modelling: towards optimal 
interpolation 



o 

(N 

> 

Q 

o 



Oh; 

o 



> 

in 



Joel Berge^'^'^*, Sedona Price^f , Adam Amara^ and Jason Rhodes^'^ 

^Department of Physics, ETH Zurich, Wolfgang- Pauli-Strasse 27, CH-8093 Zurich, Switzerland 

^ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, MS 169-327, Pasadena, CA 91109, USA 
^California Institute of Technology, 1200 East California Blvd, Pasadena, CA 91125, USA 



Accepted ???. Received ???; in original form ??? 



ABSTRACT 

Point Spread Function (PSF) modeling is a central part of any astronomy data anal- 
ysis relying on measuring the shapes of objects. It is especially crucial for weak grav- 
itational lensing, in order to beat down systematics and allow one to reach the full 
potential of weak lensing in measuring dark energy. A PSF modeling pipeline is made 
of two main steps: the first one is to assess its shape on stars, and the second is to 
interpolate it at any desired position (usually galaxies). We focus on the second part, 
and compare different interpolation schemes, including polynomial interpolation, ra- 
dial basis functions, Delaunay triangulation and Kriging. For that purpose, we develop 
simulations of PSF fields, in which stars are built from a set of basis functions defined 
from a Principal Components Analysis of a real ground-based image. We find that 
Kriging gives the most reliable interpolation, significantly better than the tradition- 
ally used polynomial interpolation. We also note that although a Kriging interpolation 
on individual images is enough to control systematics at the level necessary for current 
weak lensing surveys, more elaborate techniques will have to be developed to reach 
future ambitious surveys' requirements. 

Key words: gravitational lensing - methods: data analysis - methods: statistical 



O 1 INTRODUCTION 

I ' An accurate model of the Point Spread Function (PSF) 
is a crucial step in astronomical analyses relying on the 
estimation of galaxy shapes. For instance, an imperfect 
. !^ " PSF model has been identified as a prominent system- 
I atic effect in weak gravitational lensing measurements (e.g. 
■ iHevmans et al.ll200g : lMassev et al.ll2007l ), an d much effort is 
Ci ' unde r way to optimal l y mod el th e PSF (e.g . iKitching et al.l 
I2OIOI ). IJarvis fc Jaiiil |2004|) and |jee et al.i l|200'<t ) have de^ 
veloped techniques based on Principal Components Anal- 
ysis (PCA) to describe the PSF shape. Using an analyt- 
ical approach and an empirical one based on shapelet s 
iMassev fc RefregieJ l2005l ). IPaulin-Henriksson et all (|2008l ') 
have investigated the impact of imperfect shape descrip- 
tion of the PSF on dark energy constraints, and set re- 
quire ments on the number of stars needed to calibrate the 
PSF. lPaulin-Henriksson et al.l (|2009l ') explored how the spar- 
sity and the complexity of a PSF model affect the level of 
systematics in weak lensing surveys. Recognizing that any 
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modeling method -relying on the interpolation of the incom- 
plete information contained by individual stars, each one 
being an imperfect realization of the PSF- can be reliable 
only on scales large r than the typical distance between stars, 
lAmara et al.l l|2010l ') investigated the hardware-software bal- 
ance: although the PSF can be corrected for on large scales 
through thorough modeling, it is hardware driven on small 
scale, requiring that the telescope's characteristics are set 
such that th e PSF is kno wn, stable and well controlled at 
small scales. iRowd l|2010l l has discussed how to diagnose a 
given PSF model. In particular, he showed how measuring 
the correlation function of the residuals' ellipticity is an ul- 
timate test to the PSF interpolation. 



A first step in PSF modeling is to characterize it where 
it is sampled (i.e. on stars). Several quantities can be used 
to describe the PSF properties, such a s its moments (e.g. 
iKaiser et all 1 19951 : [Rhodes et all |2000|). or its full shape 
through shapelets coefficients ( Berge et al.i ,2008) or Prin- 
cipal Components (|jee et al.ll2007l '). [jee et all (120071 ') have 
investigated how wavelets, shapelets, and Principal Com- 
ponents perform to measure the stars' shapes. They found 
that among those three techniques. Principal Components 
provide the best description of the stars' shapes. We thus 



© ??? RAS 



2 Joel Berge et al 



follow their method to extract the shape information from 
stars. 

The next step is to interpolate the stars' shape informa- 
tion to the positions of interest, usually galaxies. The goal of 
this paper is to assess which interpolation scheme performs 
best to model a realistic PSF field, such as what can be 
dealt with in weak lensing analyses. To this end, we develop 
simulations of PSF fields, where stars are built from a set of 
basis functions defined as the Principal Components of real 
stars in a Subaru image. We then run our PSF modeling 
pipeline, including a PCA and an interpolation of the stars' 
principal components coefficients, with different interpola- 
tion schemes. 

We shall review the different interpolation schemes that 
we use in Sect. [5] Section [3] presents the simulations that we 
create and use to compare the behavior of those interpola- 
tion schemes. Our results are shown in section [l] We discuss 
them, including the dependence of the interpolation schemes 
on the stellar density, their sensitivity to outliers, and their 
impact on weak lensing systematics, in Sect. (5] We conclude 
in Sect. El 



2 INTERPOLATION SCHEMES 

This section presents the basics of the interpolation schemes 
compared in this p aper. More details can be found in the 
Numerical Recipes l|Press et all 1 19921 ). iPlaten et all (120111 ) 
proceed with similar comparisons, in the context of nonlin- 
ear density field reconstruction; they give additional useful 
mathematical background to some techniques used here. 



2.1 Polynomial interpolation 

It is common practice in weak lensing to assume that 
the PSF varies smoothly across the field, and to inter- 
polate it from st ars wit h a bivariate polynomial (e.g. 
Van Waerbeke et a l. 2005; Mivazaki et al.ll2007l : lBerge et al.1 
20081 : iFu et all |2008| ). Different variations have been in- 
vestigated and used on real data: rational fractions 
l|Van Waerbeke et al.ir2005l ). decomposition of the PSF field 
into subfields corresponding to each chip of the camera to 
help capture each chip's intrinsic behavior, or superposition 
of bivariate polynomials (Barney Rowe, private communica- 
tion). 

At the core of polynomial interpolation lies the assump- 
tion that the PSF spatial variation can be described by a 
smooth analytic function (for instance, a polynomial). How- 
ever, there is no strong physical reason why the PSF field 
should vary as a polynomial (e.g., the intrinsic camera's PSF 
can be discontinuous from chip to chip). Hence, a simple bi- 
variate polynomial is frequently no more than a rather good 
approximation of the PSF field, which may not capture de- 
tails well enough for precision-cosmology. Moreover, one usu- 
ally does not need to know such an analytic form, but only 
needs to interpolate the PSF characteristics at positions of 
interest (e.g. galaxies). It is therefore completely legitimate 
to consider the PSF field as a spatial random field, and to 
turn to less- or non-analytic forms of interpolation, such as 
those presented below. 



2.2 Radial Basis Function interpolation 

The Radial Basis Function (RBF) interpolation assumes 
that the spatial variations of the field y(x) can be repre- 
sented by the superposition of local functions 4'{r), which 
depend only on the distance to data points j, r = |x — Xj|: 



(1) 



JV-l 

2/(x) = ^tJ;i(?!)(|x-Xi|) 

where A*' is the number of data points and Wi are un- 
known weights, estimated from the value of the fields on 
data points. 

2.3 Shepard interpolation 

The Shepard interpolation is a special case of a RBF inter- 
polation where (j^ir) — !> oo as r — !> 0. It can be shown that in 
this case: 



Zi^o y(xi)0(|x-Xi| 



(2) 



In the following, we test a Shepard interpolation with 0(r) = 
r^*", with p> Q. 

2.4 Kriging interpolation 

Kriging is a Gaussian process regression technique which has 
been shown to pr ovide the best linear unbiased estimator of 
a statistical field (|Cressielll98"i . Il993l ). 

We assume that we have A'^ data points Xi, where the 
field Hi — y(xi) is known, and that we want to estimate the 
field at a given point xq. Kriging looks for the value of the 
field at this position as a weighted linear combination of the 
nearby values at known positions: 



JV-l 



(3) 



The weights Xi are estimated such that they minimize 
the error with respect to the data according to the mean 
square variation. Therefore, Kriging relies on the estimation 
of the variogram of the data to interpolate, the variogram 
being the mean square variation of the values of the field 
j/(x) as a function of the offset distance r: 



(r)~i([y(x + r)-y(x)]^) 



(4) 



where the average is over all x and r. In the following, under 
the assumption of isotropy, we assume that the variogram 
only depends on the distance r = |r|. 

We note Vij = ^dxi — Xj|) and VQj ~ v{\^g — Xj|). If 



we define the vectors 

Y = (yo,yi, . • . ,yjv-i,0) and 

Vg = {VG,0, VG,!, ■ ■ ■ , VG,N-1, 1) 

as well as the matrix 



/ voo 



V 



Vll 



«JV-1,0 Vn-1,1 

1 1 



V0,N-1 i 
Vl,N-l 1 

'^JV-l.JV-l 1 



(5) 
(6) 



(7) 



0/ 
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then the Kriging interpolation estimate yc ~ J/(xg) is 

ya = VgV-'Y. (8) 

The extra row and column in V, as well as the extra 
elements in Y and Vq allow the estimator to be unbiased. 

We find that an exponential variogram, defined as 
v{r) = exp{—r /a), where is the data points' variance 
and a is a free parameter, gives very competitive results. The 
free parameter a is the range, which sets how far from the 
position of interests data points should be used; at distances 
larger than a, the variogram becomes constant. Another free 
parameter can be added, called the nugget, which describes 
the degree of correlation of the field at very small scales, 
vo — limr-+o«(r); in this work, we find that setting t;o = 
(i.e., assuming the field does not decorrelate at very small 
scales) gives the best results. 

2.5 Delaunay triangulation 

Delaunay triangulation is a particular type of triangulation. 
Drawing triangles whose vertices are the data points, the 
Delaunay triangulation is that which maximizes the trian- 
gles' angles and minimizes their sides' length. The function 
to be interpolated is known at the triangles' vertices, then 
interpolated in their interior, where it must be estimated. In 
this paper, we use a linear interpolation between the trian- 
gles' vertices. 



3 SIMULATIONS 

In order to test the different interpolation schemes presented 
in section[21 we developed versatile simulations. They create 
mock PSF fields based on real images, whose spatial vari- 
ation, star density and star signal-to-noise ratio (S/N) can 
be set. 

Our goal is to create mock PSF fields that are as real- 
istic as possible, and can be tailored to images from current 
surveys. We thus base our simulations on real images, that 
we call 'reference images'. To create a simulated image, we 
first ext ract the stars from the r eference image using SEx- 
tractor (|Bertin fc Arnout3ll996l l. We define a set of basis 
functions, with which a mock star can be created, by pro- 
ceeding to a Principal Components Analysis of the stars 
from the real image; the associated Principal Components 
(PC) eigenvectors define the set of basis functions. During 
this process, we keep 90% of the variance information of the 
stars, so that the shape information is well extracted, while 
the background is discarded. 

We then make simulated stars by adding linearly the 
weighted PC eigenvectors. To obtain various shapes of stars, 
the weights are chosen randomly, with the only constraint 
that the variance of any given weight is similar to that of the 
coefficient associated with the corresponding PC eigenvector 
of the real stars, to ensure that mock stars are realistic. 
Finally, we require that the PSF varies across the image. 
We set its spatial variation by defining the spatial variation 
of each weight of the PC reconstruction: the random field 
from which we draw the weights is therefore given a re alistic 
power spectrum, following the technique described bv lRowd 
l|2010h . 

For the sampling of the PSF field, i.e. the position of 



Kstars stars across the image, we can either (1) use the po- 
sition of stars of the real image (in this case, ristars is the 
number of stars in the real image) or (2) lay any number 
Wstars of stars at random positions. The former possibility is 
especially useful to diagnose the PSF model that we can ex- 
pect to obtain when analyzing the real image; for instance, 
any sub-populated area may suffer from a badly constrained 
model, which will be revealed when analyzing the mock. The 
latter is useful for general modeling method development, 
e.g. to test the behavior of a given method on the stellar 
density. We use this option in this paper. 

Besides stars, we simulate rigais extra PSFs at random 
positions, that we call 'galaxy-PSFs' in the remainder of 
this paper. They are at the positions at which we want to 
interpolate the PSF from the ristars above. Simulating those 
PSFs, instead of only storing their characteristics, allows 
us to have a realistic image, on which we can run our PSF 
measurement pipeline in the same way as we do on real data 
(that is, there is no need to tailor it to the simulations), as 
well as to directly compare the PSF image and the model, 
and to use the exact same techniques to measure the PSF's 
characteristics on the simulated image and on the modeled 
image. 

A background with the same statistics as that of the 
reference image is then added to the simulated image. The 
S/N distribution of simulated stars is set by hand, either 
following that of the real image, or independent of it, to test 
the behavior of the interpolation technique on stars' S/N. 

In this paper, we base our analyses on 50 mock Sub- 
aru images, with a reference image being a typical Sub- 
aru SuprimeCam image. The simulations are 10,032x7,769 
pixels'^ in area, corresponding to 33x26 arcmin^. The PC 
decomposition of the stars in the reference image provides 
us with a set of 13 basis functions, which are used to make 
simulated stars. We set the spatial variation of the simu- 
lated stars' PC weights with a power spectrum of slope 11/4, 
which is close to what is measured for ground-based PSFs 
from the Kolmogorov spectrum for atmospheric turbulence 
(e.g. Sasicla 1994). The resulting PSF's ellipticity is of av- 
erage 0.05, and of standard deviation 0.02, consistent for all 
simulations. Each simulation contains 1000 stars (where the 
PSF would be sampled for a real image's analysis) randomly 
distributed, all with the same S/N=100. This stellar density 
of 1.1 star per arcmin^ is close to what is usually measured in 
real surveys. 'Galaxy-PSFs', where the PSF must be inter- 
polated, are added in two versions of each simulation: one 
version features randomly positioned 'galaxy-PSFs', while 
'galaxy-PSFs' of the other version are on a regular grid. The 
latter version allows us to use powerful 2D plots; the former 
is closer to reality and is used to measure correlation func- 
tions. We input 2000 'galaxy-PSFs' in each version of each 
simulation. In the following, only those 'galaxy-PSFs' are 
used to assess the correctness of any interpolation. Figure [T] 
shows two PSFs in one of our simulations. The left panels 
of Fig. [2] show a PSF ellipticity field on the 'galaxy-PSFs' 
grid (top) and on 'stars' (bottom). 

3.1 Quantities of interest 

We consider two quantities of interest to characterize the 
simulated PSF and test the model: Principal Components 
coefficients and Gaussian-weighted ellipticities 



© ??? RAS, MNRAS OOO.ITlfTSl 



4 Joel Berge et al 



■ 


, ■ 

■ 


1 


1 




■ 




1 



■T — " 



.1 . 



.■1 . -I ^ ■ ■ •■I ■■" ■ 1 




Figure 1. Example of PSFs in one simulation, with a logarithmic 
color scale. The pixel scale is 0.2 arcsec/pixel. 



^ ei + ie2 



>ii - Q22 + 2iQi2 



(9) 



Qii + Q22 

where Qij, = 1,2) are the PSF's quadrupoles. 

To compare the modeled PC coefficients with the in- 
put ones, we force our PCA-PSF measurement pipeline to 
use the exact same set of PC as those used for the simula- 
tion. This downgrades the capacity of our pipehne, since it 
is not free to find the PC eigenvectors' set that optimaUy 
describes the PSF. However for high enough S/N stars, this 
has a negligible effect. Moreover, when defining mock stars 
from the reference image's PC eigenvectors, we are careful 
to normalize the stars in the same way as we do in our PCA- 
PSF measurement pipeline, to ensure that PC keep the same 
meaning during the entire process from simulation to mod- 
eling of the PSF to test of the model. Therefore, it makes 
sense to compare the interpolated PC coefficients with those 
input in the simulations. 



4 RESULTS 

We run our PSF modeling pipeline on the simulations pre- 
sented in the previous section. As mentioned above, we force 
our PCA analysis to use the same PC eigenvectors as those 
used to create the simulated stars. Doing so, our PCA de- 
composition of an infinite S/N star would output the exact 
input star. Hence, two sources of error are present when we 
compare the interpolated PSFs with those input in the sim- 
ulations: the error coming from the shape measurement of 
finite S /N stars, and that of the interpolation. In this paper, 
we are concerned with estimating the latter. Having high 
enough S/N stars, we checked that the errors from the shape 
measurement are small enough that they can be neglected 
in our first tests (Fig. I2I6[) . The last test, using correlation 
functions (sect. [4^ is sensitive to the errors from the shape 
measurement; we will use them to discuss the errors coming 
from the interpolation. 

Our results are presented below. 



Polynomial interpolation We interpolate the PSF's PC 
coefficients with: (1) a bivariate polynomial on the entire 
field, the degree of which is focused with a minimiza- 
tion, but is constrained to remain small enough to avoid fast 
divergences; and (2) 2-dimensional Chebyshev polynomials 
interpolation. 

We find Chebyshev polynomials, although they are 
bound and should be expected to give better results, not to 
perform better than regular bivariate polynomials. There- 
fore, in the remainder of this paper, we will only report 
results from regular bivariate polynomials. We find the best 
fitting polynomial to be of 5th order. 

RBF interpolation We try Gaussian-RBF, multiquadra- 
tic-RBF, inverse multiquadratic-RBF, and thin-plate-RBF. 
Although Gaussian-RBF behave better than the others, we 
find all those kinds of RBF to be too unstable, and need 
too fine a tuning to give reliable models. Therefore, we do 
not spend more time on RBF interpolations here. However, 
we find the Shepard interpolation to perform as well as a 
polynomial interpolation, with almost no fine tuning needed, 
for p ^ 3. 

Kriging interpolation We try exponential and spherical 
variograms, with varying range and nugget. We find the ex- 
ponential variogram to be more stable and to perform better 
than the spherical variogram. Therefore, in the remainder of 
this paper, we will only consider an exponential variogram 
when using Kriging interpolation. 

We get better results when setting the nugget to 0. 
When normalizing the images' coordinates so that < 2; < 1 
and < y < 1, we find that the interpolation gives better re- 
sults when setting the range of the variogram between 0.25 
and 0.75. The range acts on the stability of the model at 
small scales by taking more or less stars into account to 
estimate the variogram. 

In the following, we compare the results of the interpo- 
lation with a polynomial, a Delaunay triangulation, and a 
Kriging technique. 

4.1 Recovery of PSF's field information 

We start by comparing the Gaussian-weighted ellipticity of 
the interpolated PSFs with that of the simulated PSF, for 
a given simulation. Although the limited space prevents us 
from showing the results for several simulations, we veri- 
fied that our results are consistent for all our simulations. 
The top-left panel of Fig. [2j as seen in landscape, shows 
the simulation's ellipticity field, sampled at the position of 
'galaxy-PSFs'. The bottom-left panel show the PSF's ellip- 
ticity of the stars used for the interpolation. The middle and 
right columns of the figure show the models and residuals 
we obtained with a 5th-order-polynomial (top), Delaunay 
(middle) and Kriging (bottom) interpolations, on the grid 
formed by 'galaxy-PSFs'. 

Figure |3] shows the same information as Fig. [51 decom- 
posed into the two components of the ellipticity, ei (top) and 
62 (bottom). For both components, the top-left panel shows 
the input, and the right two columns show the ellipticity of 
the interpolated PSF, as well as the residuals, for the three 
interpolation techniques we compare here. 
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Figure 2. Unweighted ellipticity of the simulated PSF field. Top-left: ellipticity of the input PSF on a regular grid. Bottom-left: ellipticity 
of the input PSF at the position of 'stars' used as samples of the PSF for the interpolation of its Principal Components. Middle column: 
recovered PSF ellipticity fields after interpolation of the PSF PCs, for polynomial interpolation (top), Delaunay triangulation (center) 
and Kriging interpolation (bottom). Right column: residuals, amplified by a factor 3 for better visibility. 



© ??? RAS, MNRAS OOO.ITlfTsI 



6 Joel Berge et al 





2000 4000 6000 8000 2000 4000 6000 8000 



Figure 3. Same as Fig. [2] for tiie two components of the unweighted ellipticity of the PSF. Top: ei. Bottom: 62. Due to gridding aspects, 
it was not possible to represent the ellipticity components for individual stars. Note that the color scale for the residuals is narrower than 
that of the input and of the model (-0.03 to 0.03 instead of -0.1 to 0.1). 



© ??? RAS, MNRAS OOO.nifTSi 



PSF interpolation 




2000 4000 6000 8000 




2000 4000 6000 8000 2000 4000 6000 8000 




2000 4000 6000 8000 2000 4000 6000 8000 



Figure 4. Same as Fig. [2] for the first two Principal Components coefficients of the PSF. Note that the color scale for the residuals 
narrower than that of the input and of the model (-0.006 to 0.006 instead of -0.03 to 0.03). 



© ??? RAS, MNRAS OOO.ITlfTsI 



8 Joel Berge et al 





Figure 4. (Ct'd) Third and fourth Principal Components coefficients of the PSF. 
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Figure [l] with the same layout as Fig. [3l goes one step 
further in decomposing the information. It shows the inter- 
polation of the first four Principal Components (the most 
significant) of the PSF themselves. Hence, while the previ- 
ous figures showed the combined 'physical' information of 
the ellipticity, this figure is at the core of our interpola- 
tion pipeline, since it shows how each individual interpo- 
lation performs. Similar figures, made from regularly grid- 
ded 'galaxy-PSFs' simulations, cannot be done for real data, 
where one knows the PSF information only on stars. On real 
data, one can only rely on testing the interpolation of the 
ellipticity on stars (i.e. one should transform Fig. [5]such that 
information appears at the position of stars only). 

A bivariate polynomial interpolation fails to capture the 
quickly varying features of the PSF pattern: significant resid- 
uals are seen by eye around the most significant gradients 
of the PC and ellipticity fields, showing how this interpo- 
lation mostly smoothes the information. Higher and higher 
order polynomials should manage to capture increasingly 
small variations of the field; however, this requires an un- 
comfortable tuning of the polynomial, which is hardly com- 
patible with the automatization that is necessary to deal 
with big surveys. The Delaunay triangulation, although suf- 
fering non- negligible residuals when the field varies quickly 
(see e.g. the third PC), produces satisfactory interpolations. 
Finally, the Kriging technique recovers the information al- 
most perfectly, with homogeneous, near-zero residuals. Only 
when the gradients of the fields are very strong does it leave 
an imprint in the residuals maps. 

From these three figures, we can already conclude that 
the Kriging interpolation performs best, slightly better than 
a Delaunay triangulation interpolation, and significantly 
better than a bivariate polynomial interpolation. Moreover, 
the Delaunay triangulation and Kriging interpolations need 
only minimal tuning, and hence are more portable and 
adapted to diverse data sets. 

In Fig. [5l we compare the distribution of the ellipticity 
residuals, between the three considered interpolation tech- 
niques, for 50 simulations. Figure [S] provides similar com- 
parisons for the first four Principal Components of the PSF 
shape. These figures, showing that the distribution of resid- 
uals is clearly tightest for the Kriging interpolation, confirm 
our claims above about how much better this technique per- 
forms. Furthermore, no significant correlation is visible be- 
tween the residuals obtained with a polynomial interpolation 
and the other two interpolation schemes; they perform well 
even where a polynomial interpolation performs badly. On 
the other hand, the residuals obtained with a Kriging and 
those obtained with a Delaunay triangulation interpolations 
are correlated. This shows that since they are both sensitive 
to small scale information, their behaviors are comparable. 
For instance, they suffer similar limitations when the field 
varies quickly (as shown in Figs. [2]|4]). 

4.2 Ellipticity correlation functions 

We now turn to the PSF ellipticity correlation functions. If 
the interpolation is correct, the correlation function of the 
residuals should be consistent with zero, or at least negli- 
gible compare d with that of the measured PSF ellipticity. 
In particular, iRowd (|2010h has shown that the correlation 
function of the residuals helps diagnose an under- or an over- 
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Figure 5. 2D distributions of the ellipticity's residuals in the 
three planes defined by the one-to-one comparisons of the three 
interpolation schemes. Upper panel: ei. Lower panel: 62. Contours 
show the distribution of the residuals, normalized such that the 
maximum is 1. Contours start from 0.1 and increase by steps of 
0.1. 



fitting, and therefore is recognized a perfect test to choose 
between different interpolations. 

The ellipticity correlations functions are defined as 
Xd^l) e(r)e*(r + S) >, where e is the complex ellip- 
ticity of either the input simulation, the interpolated PSF, 
or the residuals, and e* is the ellipticity's complex con- 
jugate. In case of the residuals, the correla tion function 
is identical to the function Di (9) defined by I Rowel (|2010l ) 
{Di{8) =< 5e{r)*5e{r + 6) >, where 5e are the residuals). 
Although on real data, the PSF ellipticity correlation func- 
tions can only be estimated with stars, here we measure 
them using only the 'galaxy-PSFs' in our simulations, that 
is, the measurement of the PSF on points not used to per- 
form the interpolation. Hence, we discard the stars, which 
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Figure 6. 2D distributions of the first four PC coefficients' residuals in tlie three planes defined by the one-to-one comparisons of 
the three interpolation schemes. From left to right and top to bottom: first, second, third and fourth coefficient. Contours show the 
distribution of the residuals, normalized such that the maximum is f . Contours start from 0.1 and increase by steps of 0.1. 



prevents our estimation of the residuals' correlation function 
from being underestimated due to a perfect interpolation on 
stars. Our measurements of the correlation functions there- 
fore only probe the reconstruction of the PSF where it is not 
known a priori (and used for the interpolation). 

Figure [T] shows the PSF ellipticity correlation functions, 
averaged over 50 simulations, for a bivariate-polynomial, a 
Delaunay triangulation, and a Kriging interpolations from 
left to right. The black solid lines represent the ellipticity 
correlation functions as measured on the simulations. The 
red lines are those of the interpolated PSF. The green lines 
are those of the residuals. The correlation function of the 
residuals in Fig.[7|confirm that a Kriging interpolation yields 
significantly better results than a polynomial interpolation. 
The difference is less striking with a Delaunay triangula- 
tion, although clearly noticeable. We discuss the effects on 



systematics in the context of weak gravitational lensing in 
section [5TT] 



5 DISCUSSION 

In this section, after discussing the level of systematics from 
the PSF interpolation for weak lensing, we discuss the sta- 
bility of the techniques presented above. In particular, we 
focus on their dependence on the stellar S/N and stellar 
density. 

5.1 Weak gravitational lensing systematics 

Although considered a premier cosmological probe, weak 
lensing is such a small effect that all possible systematic ef- 
fects must be tackled very thoroughly. The PSF is the most 
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Figure 7. PSF ellipticity correlation functions, averaged over 50 simulations. Left: polynomial interpolation. Center: Delaunay triangu- 
lation interpolation. Right: Kriging interpolation. In each panel, the black line represents the ellipticity correlation function of the input, 
the red line shows that of the model, and the green one shows that of the residuals. Dotted lines show the upper limits on the residuals 
that one must satisfy so that the weak lensing systematics crj from the PSF are less than 10~®, 10~® and 10~^ from top to bottom. 



significant systematic, malcing the correction of shapes for 
the PSF eflects of paramount importance. Here, we quantify 
the systematics on the cosmic shear power spectrum due to 
PSF interpolation errors using the three schemes mentioned 
above. 

We assume that the PSF deconvolution, which oc- 
curs after its interpolation when measuring galaxy shapes, 
does not bring additional systematics, so that all sys- 
tematics attributable to the PSF come from its model- 
ing only. We compute the impact on the shear correl ation 
function by following the defi nition by iRowd l|2010h and 
IPauhn-Henriksson et all lj2009( ): 



mm < 



(10) 



where Di is the correlation function of the residuals (see 
above), is the shear susceptibility, and J?psp and J?gai 
are the sizes of the PSF and of ga l axies. We take the 
same values as those used by iRowd l|2010l ): P'' = 1.84, 
< (i?PSF/-Rgai)'' >~ (1/ 1.5)"^. Note that in Eg. lITOl). w e 



ignore the second term of IPaulin-Henriksson et ahl 1 2009[ )'s 
Eq. (15), which depends on the amplitude of the pre- 
correction PSF ellipticity. Indeed, since it puts the emphasis 
on the pre-correction PSF, this strong dependence makes it 
less relevant to our discussion. 

The impact on the shear correlation function (5^[jl can 
be integrated in Fourier space, to estimate the lev el of sys- 
tematics, in the sense of lAmara fc Refregiej (|2008l ). coming 
from the imperfect PSF interpolation: 



^_ / |crWf + l)dlnl 



(11) 



The dotted lines in Fig. [7] show the upper limits on 
the residuals ellipticity correlation function (assumed to be 
a constant fraction of the shear correlation function on all 
scales) such that the systematics due to the imperfect PSF 
modeling are less than aly^ = 10~^, 10~®, 10~^ from top to 
bottom. We use a ACDM cosmology with erg = 0.8 to esti- 
mate those limits. Due to our neglecting the dependence of 
CTgys on the pre-correction PSF ellipticity, the limits plotted 



in Fig. [7] are optimistic compared to those we could achieve 
in a real survey, but give a good general sense of the required 

leve l of the PSF residua ls. 

lAmara fc Refregieij (|2008l ) investigated the allowed val- 
ues of cTgys such that, for a given survey, errors in the esti- 
mation of cosmological parameters are dominated by sta- 
tistical errors. Table [1] gives the maximum crfyg allowed 
for some current a nd p lanned surveys, computed using 
lAmara fc Refreeierl (j2008l Vs scaling relation (21). By com- 
paring this table with Fig. [71 while it appears that a sim- 
ple bivariate polynomial interpolation does not allow one 
to reach the limit on the PSF systematics on current sur- 
veys such as the CFHTLS Wide, a Kriging approach lowers 
the residuals enough to meet the requirements. However, the 
techniques explored in this paper would not allow us to lower 
the systematics to the level needed for future ambitious wide 
field surveys. 

Nevertheless, our analysis is made from mock ground- 
based images, and our techniques are restricted to single- 
field interpolation. More complex PSF interpolation schemes 
can be thought of. For instance, a multi-field interpola- 
tion scheme would allow one to look for coherent pattern 
from image to image, and therefore improve on the PSF 
model by taking into account more information than avail- 
able in a si ngle image. Exampl es of such tech niques are 
provided by Ijarvis fc JainI l|2004 ) and |jee fc TV son (20111); 
a similar approach has been undertaken in COSMOS by 
ISchrabback et all (|2010l l. To sum up, our conclusions stand 
for PSF modeling based on single-field interpolation of the 
PSF from ground-based imaging surveys. More elaborate, 
multi-field interpolation techniques, either on ground-based 
or on space-based data, should allow one to lower the level 
of residuals, and reach the requirements for Stage IV weak 
lensing surveys. 



5.2 Impact of stars' S/N — shape measurement 
errors vs interpolation errors 

Although so far we have only focussed on the interpolation 
of the PCA coefficients, our analysis includes the first step 
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Table 1 . Maximum level of systematics on shear power spectrum 
for weak lensing surveys 



Survey 


Area 
(deg2) 


Galaxy density 
(arcmin"'^) 


Median 
rcdsliift 


""sys 


CFHTLS 


170 


12 


0.7 


2 X 10"'^ 


COSMOS 


2 


50 


0.9 


8 X lO"'^ 


Stage IV 


15000 


30 


0.9 


1.2 X 10-'^ 



10-^ 




0.1 1.0 10.0 

8 (arcmin) 



Figure 8. EUipticity correlation functions, with a Kriging in- 
terpolation, for different S/N stars. The black and red lines are 
the input and model correlation functions for S/N=100 stars. The 
green line is the correlation function of the residuals for S/N=100 
stars. The blue line is the correlation function of the residuals for 
S/N=700 stars. The shaded region in the upper part is the re- 
gion spanned by the correlation functions of theoretical Gaussian 
random fields used to estimate the optimal correlation function 
of residuals after Wiener-filtering (lower shaded region). 

of a PSF modeling process, namely the PSF's shape mea- 
surement (PC A decomposition in our pipeline). 

The non-vanishing ellipticity correlation functions of 
the residuals (Fig. (Tjl betray a dependence of the PSF model 
on the stellar S/N, defined as the ratio of the mean of the 
star's pixels to the r.m.s of the background. To test this ef- 
fect of the S/N, we create a similar set of simulations, with 
higher S/N stars (S/N=700 for each star). The blue line in 
Fig. E] shows the correlation function of the residuals for a 
Kriging interpolation with those higher S/N stars, and com- 
pares it with the residuals obtained previously with stars of 
S/N=100 (green line). It is obvious that the residuals have 
been shifted downwards, due to the better shape measure- 
ment in the first step of our pipeline. 

To discriminate the errors from the shape measurements 
against those from the interpolation, we estimate the min- 
imal errors on the PSF modeling that we can expect in 
the ideal case of a Wiener-filtered Gaussian random field 
ijAmara et al In this case, the model is optimal (i.e., 

it gives the lowest achievable residuals), and depends only 



on the number of stars and on their S/N. The shaded region 
in the upper part of Fig. [8] shows the span in ellipticity cor- 
relation functions of 12 Gaussian random fields with power 
spectrum resembling that of our simulations. The shaded 
region in the lower part of the figure is the span of the cor- 
responding correlation functions of the residuals, for stars 
with S/N=100. At small scales, the correlation function of 
the residuals reflects the errors from the shape measurement. 
For increasing scales, the interpolation cancels out those er- 
rors, and the residuals tend to zero at large scales. Although 
these perfect theoretical expectations cannot be rigorously 
compared with our simulations, since they are made in a 
different regime (Gaussian vs non-Gaussian flelds interpo- 
lated with different techniques), they give a good indica- 
tion about how the errors from the interpolation compare 
with those from the shape measurement in our analysis. At 
small scales, the order of magnitude of the correlation func- 
tion of the residuals (green line) being similar to that of 
the theoretical estimation, we can safely conclude that our 
Kriging interpolation is likely close to optimal. Our theoret- 
ical estimates do not allow us to explain the presence of a 
plateau in the correlation function of the residuals at large 
scale, since at large scales, Kriging is underperforming com- 
pared with theoretical expectations. But the increase of the 
pre-correction ellipticity correlation function at large scales 
in our simulations, which is not taken into account in our 
theoretical estimations, may explain the presence of these 
plateaus. 

Finally, we tested the robustness of the interpolation 
techniques against the stellar S/N (and hence, the reliability 
of the stars' model) by creating another set of simulations, 
with a realistic S/N distribution, ranging from S/N=20 to 
S/N=1000, and peaking around S/N=100. The total num- 
ber of stars is unchanged. Lowering the S/N of some stars 
decreases the fidelity of the associated PCA model, and po- 
tentially perturbs the interpolation, since the coefficients on 
which it is based are more noisy. We find that ignoring low 
S /N stars does not impact the interpolation, as long as high 
enough S/N (^ 100) stars are suf ficiently numerous. This 
result was already mentioned by IPaulin-Henriksson et al.l 
(|2009l ) as part of their work on PSF shape measurement. 
They claimed that ~ 1 star per square arcminute with 
S/N^lOO was necessary to obtain a reliable model. 

5.3 Global vs local stellar density 

Here, we investigate the dependence of the interpolation 
techniques on the number (identically, the density) of stars. 
We compare the fidelity of the interpolation at a given po- 
sition, as a function of the local stellar density and as a 
function of the global stellar density in the image. In the 
former case, we naively expect interpolations to be better 
constrained in high density regions than in lower density re- 
gions. In the latter case, we expect interpolations to perform 
better with more stars. 

Figure|9]shows box-plots of the distribution of the resid- 
uals of the two ellipticity components, as a function of the 
local star density, computed in one-arcmin-radius circular 
apertures around 'galaxy-PSFs', for a polynomial interpola- 
tion (black), a Delaunay triangulation (red) and a Kriging 
interpolation (green). Each distribution is drawn from the 
set of 50 simulations used in section |4] For each box-plot. 



© ??? RAS, MNRAS OOO.ITlfTSl 



PSF interpolation 13 



the central line represents the median of the distribution, the 
boxes around it are the second and third quartiles, and the 
external symbols are the outer 5% points. Figure [TOl shows 
the r.m.s of those distributions, as a function of local stellar 
density. 

These figures concur with our claims above, that the 
Kriging interpolation gives the smallest residuals. This is 
the case for all stellar densities. They also show that the 
variance of the residuals, as well as central quartiles of their 
distributions, does not exhibit a notable dependence on the 
local stellar density. Therefore, high density regions do not 
constrain the interpolation more than lower density regions. 
However, it is clear from Fig. |9] that the width of the full 
distribution of residuals decreases in high stellar density re- 
gions: outliers from interpolation failures are much less likely 
in high stellar density region than in low density regions. 
This observation can provide a way to weigh the PSF model 
in an image. 

Figure [11] shows a similar box-plot as Fig. |9l but for 
ellipticity residuals as a function of the global stellar den- 
sity (i.e., as a function of the total number of stars available 
in the image to perform the interpolation). The accuracy of 
a polynomial interpolation does not show a strong depen- 
dence on the total number of star. This is because, for a low 
enough order polynomial, the general pattern of the poly- 
nomial is set by a small number of stars; adding stars does 
not help constrain the polynomial better. On the opposite, 
a Delaunay triangulation and Kriging use all information at 
small scales: their dependence on the total number of stars 
is significant. 

To sum up this section, the interpolation accuracy 
mostly depends on the global star density rather than on the 
local star density, especially for Delaunay triangulation and 
Kriging. While the global star density acts on the general 
accuracy of the interpolation, the local star density helps 
eliminate outliers from the interpolation. 

5.4 Sensitivity to outliers 

Each star in an image provides an imperfect realization of 
the PSF, which is pixelated and noisy. Furthermore, the PSF 
shape measurement does not give a perfect description of 
each star. From those facts, outliers are possible, which could 
be seen e.g. as stars with unphysical ellipticity. Such outliers 
can be source of error in the interpolation. For instance, one 
can see a small outlier in the stars' ellipticity on Fig. O at 
position {x,y) « (1500,5500). It is visible that the Delau- 
nay triangulation and the Kriging interpolation are affected 
by this outlier, since their residuals are clearly bigger than 
average around this position. 

We find that, due to their being local interpolation 
schemes, the Delaunay triangulation and the Kriging inter- 
polation are more prone to outliers than a polynomial in- 
terpolation, which tends to smooth the spatial pattern out 
(however, a very high order polynomial interpolation would 
also become more prone to outliers). 

This sensitivity to outliers is alleviated by a simple 
control of outliers. We perform an Interquartile range test 
on the distribution of each Principal Component. For each 
component, outliers are defined as having values less than 
Q\ — alQR or more than Q3 -I- alQR, where Qi is the ith 
quartile of the component's distribution, IQR = Q3. — Qi is 
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Figure 9. Ellipticity residuals box plots as a function of local stel- 
lar density, for polynomial interpolation (black), Delaunay trian- 
gulation (red) and Kriging (green). For each box-plot, the central 
line represents the median of the distribution of residuals, the ad- 
jacent boxes are the central quartiles, and the symbols are the 5% 
outer points. For each local stellar density, box-plots are offsets 
for visibility. 
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Figure 10. Median and rms of ellipticity residuals distributions 
as a function of stellar density. The color code is the same as in 
Fig.E] 



the interquartile range, and a is a free parameter. As soon 
as at least one component of a star is found to be an out- 
lier, this star is discarded. Using a loose criterium (a — 3.5) 
helps control the sensitivity to outliers, while removing only 
a negligible proportion of stars (typically 1%). 
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Figure 1 1 . Same as Fig. |9] but for global star density. 



6 CONCLUSION 

We investigated iiow to improve on the techniques usually 
used to interpolate the PSF. An extremely reliable PSF 
model is necessary, in particular for weak gravitational lens- 
ing analyses, to beat systematics down. Any PSF model re- 
lies on a robust interpolation of any quantity used to de- 
scribe the PSF (e.g. its moments, or the coefficients of a 
decomposition on a given set of basis functions). 

To this end, we developed simulations of PSF fields. 
The simulations are based on a typical Subaru SuprimeCam 
image, whose stars are decomposed into Principal Compo- 
nents, which are then used as basis functions to build mock 
stars. The background of the simulations is given the same 
statistics as that of the real image used as a reference for 
our simulations. 

We used the simulations to compare several interpola- 
tions schemes, including bivariate polynomial, radial basis 
functions, Delaunay triangulation and Kriging. We assumed 
that the first step of PSF modeling (its shape measurement 
-for instance, its decomposition into PCA in our analysis) 
is well constrained and brings only negligible errors on the 
model. Therefore, we focussed on the interpolation part of 
the pipeline, that is, the interpolation of the principal com- 
ponents coefficients of the stars. We found RBF to be unsta- 
ble. We found a Kriging approach to bring the best results, 
slightly better than a Delaunay triangulation, and signifi- 
cantly better than the traditionally used bivariate polyno- 
mial interpolation. 

We showed how the full PSF modeling depends upon 
the S/N of stars, and how the Kriging interpolation is close 
to optimal. We discussed how the different interpolation 
schemes depend on the local and global star densities: mod- 
els are better controlled in high stellar density regions; De- 
launay triangulation and Kriging are more sensitive than a 
polynomial interpolation to the total number of stars. We 
noticed that Kriging and Delaunay interpolations are more 
prone to outliers than a polynomial interpolation; however, 
a simple control on outliers helps cancel this difficulty. Fi- 



nally, although Kriging is numerically more expensive than 
a polynomial interpolation, it is worth the improvement it 
provides. 

The ellipticity correlation functions of the residuals be- 
tween the input PSF and the interpolated PSF allowed us 
to conclude that although a Kriging interpolation is suffi- 
ciently accurate for current weak lensing surveys, it will not 
allow us, as used in this paper, to meet the requirement on 
the control of systematics for upcoming ambitious surveys. 
Since Kriging is considered as the best linear unbiased inter- 
polator, it seems difficult to go beyond it to make significant 
improvements. Etowever, our analysis stands for single-field 
interpolation, i.e., the PSF is interpolated on each image 
separately, with a maximum density of useful stars about 
one per square arcminute. More elaborate techniques, re- 
lying on several images to perform the interpolation, will 
improve on the results of this paper. For example, when the 
PSF is know to be stable in time (i.e., for space-based ob- 
servations), stacking several images allows one to artificially 
increase the density of stars, and therefore to better con- 
strain the interpolation. Another possibility is to look for 
coherent patterns in different images, even for ground-based 
data, which allows for a multi-field interpolation to be per- 
formed. Those new techniques will surely permit to meet the 
requirements for future surveys. 
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